100 REM N-BODY SIMULATOR
110 REM VERSION 6.09
120 REM BY RICHARD LUCAS
130 :
140 REM INITIALIZE
150 POKE2038,PEEK(55):POKE2039,PEEK(56)
160 POKE 56,62: POKE 55,0: CLR
170 DIM X(40),Y(40),Z(40),U(40),V(40),W(40),X1(40),Y1(40),Z1(40)
180 DIM M(40),GM(40),E2(7)
190 DIM X0(40),Y0(40),Z0(40),U0(40),V0(40),W0(40),A0(40),B0(40),C0(40)
200 DIM EX(7),UN(3),UN$(3),TU$(3)
210 FORI=1TO3:READTU$(I),UN$(I):NEXT
220 DATA "DAYS","10^7KILOMETERS-TONS-DAYS","DAYS","[193][213]-[197]ARTH MASS-DAYS","SECONDS"
230 DATA 1000KM-KG-SEC
240 GOSUB 2820
250 G=6.67E-11:SY=1:CF=40:C4=504:C7=7:C8=248:VI=53248:HI=VI+16:C5=255
260 HR=16*1024:O$="EPNCDLS"+CHR$(133)+CHR$(136)+"XYZUVWMTR[196]":SP=0:DV=8
270 NB=0:DT=10:D2=DT*DT/2:D3=D2*DT/3:D4=D3*DT/4:CB=1
280 FORI=0TO7:EX(I)=2^(7-I):NEXT
290 FORI=0TO7:E2(I)=2^I:NEXT
300 UN(1)=1000*86400^2/1E10^3
310 UN(2)=5.9742E24*86400^2/1.4959789E11^3
320 UN(3)=1/1E6^3
330 POKE 53280,0:POKE 53281,0:PRINT CHR$(14)+CHR$(8)+CHR$(151);
340 B$="":L$="":FORI=1TO38:B$=B$+CHR$(32):L$=L$+CHR$(192):NEXTI
350 B$=CHR$(29)+B$+CHR$(29)
360 :
370 REM MAIN LOOP
380 GOSUB 2460
390 GOSUB 2620
400 PRINT ">";
410 GET A$:IFA$=""THEN 410
420 FORI=1TOLEN(O$):IFA$=MID$(O$,I,1)THEN PRINTA$:GOTO 450
430 NEXT
440 GOTO 410
450 ON I GOTO 490,510,1970,2900,3070,2110,2250,3190,3230
460 IF I=18 THEN INPUT"[206]EW NAME OF BODY";N$(CB):GOTO 370
462 IF I=19 THEN INPUT"[196]RIVE #";DV:GOTO 370
470 INPUT"NEW VALUE";NV
480 ON I-9 GOTO 2390,2400,2410,2420,2430,2440,2380,2450
490 POKE55,PEEK(2038):POKE56,PEEK(2039):CLR:END
500 :
510 REM PLOT TRAJECTORIES
520 IF NB=0 THEN PRINT "[206]O BODIES IN CURRENT SYSTEM.":GOTO 370
530 IF SP THEN 610:REM SKIP HIRES FOR SPRITES
540 REM SET UP HIRES SCREEN
550 POKE 56578,PEEK(56578)OR3:REM SWITCH TO VIC BANK 1 (16K-32K)
560 POKE 56576,(PEEK(56576)AND252)OR2
570 POKE 53272,(PEEK(53272)AND15)OR128:REM CHAR SCREEN IS IN 9TH K
580 POKE 53265,PEEK(53265)OR32:REM TURN ON HIRES SCREEN
590 POKE 820,0:POKE 821,64:POKE 822,0:POKE 823,96:POKE 251,0:SYS 49152
600 POKE 820,0:POKE 821,96:POKE 822,231:POKE 823,99:POKE 251,16:SYS 49152
610 IFSPTHENFORI=VITOHI:POKEI,.:NEXT
620 IFSPTHEN A1=0:FORI=1TO8:A1=A1 OR(-(I<=NB)*2^(I-1)):NEXT:POKE VI+21,A1
630 IFSPTHENPOKE 53281,0:PRINTCHR$(147)
640 T=0
650 :
660 REM MOVE START PARAMETERS TO WORKING ARRAYS
670 FORI=1TONB
680 X(I)=X0(I)
690 Y(I)=Y0(I)
700 Z(I)=Z0(I)
710 U(I)=U0(I)
720 V(I)=V0(I)
730 W(I)=W0(I)
740 NEXT
750 REM COMPUTE ACCEL AT TIME DT BEFORE START
760 FORI=1TONB
770 A0(I)=.:B0(I)=.:C0(I)=.
780 NEXT
790 FORI=1TONB
800 AX=.:AY=.:AZ=.
810 FORJ=1TONB
820 IFI=JTHEN910
830 DX=X(J)-X(I)
840 DY=Y(J)-Y(I)
850 DZ=Z(J)-Z(I)
860 R=SQR(DX*DX+DY*DY+DZ*DZ)
870 R3=R*R*R/GM(J)
880 AX=AX+DX/R3
890 AY=AY+DY/R3
900 AZ=AZ+DZ/R3
910 NEXT
920 X1(I)=X(I)-U(I)*DT+AX*D2
930 Y1(I)=Y(I)-V(I)*DT+AY*D2
940 Z1(I)=Z(I)-W(I)*DT+AZ*D2
950 NEXT
960 FORI=1TONB
970 AX=.:AY=.:AZ=.
980 FORJ=1TONB
990 IFI=JTHEN1080
1000 DX=X1(J)-X1(I)
1010 DY=Y1(J)-Y1(I)
1020 DZ=Z1(J)-Z1(I)
1030 R=SQR(DX*DX+DY*DY+DZ*DZ)
1040 R3=R*R*R/GM(J)
1050 AX=AX+DX/R3
1060 AY=AY+DY/R3
1070 AZ=AZ+DZ/R3
1080 NEXT
1090 A0(I)=AX
1100 B0(I)=AY
1110 C0(I)=AZ
1120 NEXT
1130 :
1140 REM CALCULATE NEW SYSTEM STATE
1150 FORI=1TONB
1160 A1=.:B1=.:C1=.:A2=.:B2=.:C2=.
1170 FORJ=1TONB
1180 IFI=JTHEN1270
1190 DX=X(J)-X(I)
1200 DY=Y(J)-Y(I)
1210 DZ=Z(J)-Z(I)
1220 R=SQR(DX*DX+DY*DY+DZ*DZ)
1230 R3=R*R*R/GM(J)
1240 A1=A1+DX/R3
1250 B1=B1+DY/R3
1260 C1=C1+DZ/R3
1270 NEXT
1280 J0=(A1-A0(I))/DT
1290 K0=(B1-B0(I))/DT
1300 L0=(C1-C0(I))/DT
1310 X2=X(I)+U(I)*DT+A1*D2+J0*D3
1320 Y2=Y(I)+V(I)*DT+B1*D2+K0*D3
1330 Z2=Z(I)+W(I)*DT+C1*D2+L0*D3
1340 FORJ=1TONB
1350 IFI=JTHEN1440
1360 DX=X(J)-X2
1370 DY=Y(J)-Y2
1380 DZ=Z(J)-Z2
1390 R=SQR(DX*DX+DY*DY+DZ*DZ)
1400 R3=R*R*R/GM(J)
1410 A2=A2+DX/R3
1420 B2=B2+DY/R3
1430 C2=C2+DZ/R3
1440 NEXT
1450 J1=(A2-A1)/DT
1460 K1=(B2-B1)/DT
1470 L1=(C2-C1)/DT
1480 M1=(A2-2*A1+A0(I))/(DT*DT)
1490 N1=(B2-2*B1+B0(I))/(DT*DT)
1500 O1=(C2-2*C1+C0(I))/(DT*DT)
1510 X1(I)=X(I)+U(I)*DT+A1*D2+J1*D3+M1*D4
1520 Y1(I)=Y(I)+V(I)*DT+B1*D2+K1*D3+N1*D4
1530 Z1(I)=Z(I)+W(I)*DT+C1*D2+L1*D3+O1*D4
1540 U1(I)=U(I)+A1*DT+J1*D2+M1*D3
1550 V1(I)=V(I)+B1*DT+K1*D2+N1*D3
1560 W1(I)=W(I)+C1*DT+L1*D2+O1*D3
1570 A0(I)=A1
1580 B0(I)=B1
1590 C0(I)=C1
1600 NEXT
1610 :
1620 FORI=1TONB
1630 X(I)=X1(I)
1640 Y(I)=Y1(I)
1650 Z(I)=Z1(I)
1660 U(I)=U1(I)
1670 V(I)=V1(I)
1680 W(I)=W1(I)
1690 X=X(I):IFX(I)<.ORX(I)>319THEN1720
1700 Y=Y(I):IFY(I)<.ORY(I)>199THEN1720
1710 GOSUB2710
1720 NEXT
1730 T=T+DT
1740 IF SP THEN PRINTCHR$(19);T
1750 GETA$:IFA$=""THEN1150
1760 :
1770 REM RESTORE CHARACTER SCREEN
1780 POKE53265,PEEK(53265)AND223
1790 POKE 56578,PEEK(56578)OR3
1800 POKE 56576,(PEEK(56576)AND252)OR3
1810 POKE 53272,(PEEK(53272)AND15)OR16
1820 POKE VI+21,0:REM TURN OFF SPRITES
1830 POKE 53281,0
1840 IF A$=CHR$(133) THEN 1860
1850 GOTO 370
1860 PRINT CHR$(17);"[211]TORING PRESENT SYSTEM IN MEMORY.":GOSUB 3160
1870 FOR I=1 TO NB
1880 X0(I)=X(I)
1890 Y0(I)=Y(I)
1900 Z0(I)=Z(I)
1910 U0(I)=U(I)
1920 V0(I)=V(I)
1930 W0(I)=W(I)
1940 NEXT I
1950 GOTO 370
1960 :
1970 REM GET NEW SYSTEM FROM USER
1980 INPUT"[206]UMBER OF BODIES";NB
1990 IFNB<1ORNB>50THEN370
2000 FOR I=1 TO NB
2010 CB=I:GOSUB 2460
2020 PRINT"[206]AME OF BODY"I;:INPUTN$(I)
2030 IFLEN(N$(I))>25THEN2020
2040 PRINT"[205]ASS OF BODY"I;
2050 INPUT M(I):GM(I)=G*M(I)*UN(SY)
2060 INPUT"[201]NPUT LOCATION IN X,Y,Z FORM";X0(I),Y0(I),Z0(I)
2070 INPUT"[201]NPUT VELOCITY IN U,V,W FORM";U0(I),V0(I),W0(I)
2080 NEXT I
2090 GOTO 370
2100 :
2110 REM LOAD SYSTEM DESCRIPTION FROM DISK
2120 PRINT "[204]OAD SYSTEM DATA FROM DISK."
2130 INPUT"[212]YPE NAME OF DATA FILE";A$
2140 IF LEN(A$)>13 THEN PRINT"[212]OO LONG.":GOTO 2130
2150 OPEN 15,DV,15:OPEN 2,DV,2,"0:"+A$+".NB,S,R"
2160 GOSUB 3170
2170 IF ER<>0 THEN PRINT ER$(1);ER$(2);ER$(3);ER$(4):GOSUB 3160:GOTO 2230
2180 INPUT#2,SY,NB
2190 FORI=1TONB
2200 INPUT#2,N$(I),X0(I),Y0(I),Z0(I),U0(I),V0(I),W0(I),M(I)
2210 GM(I)=G*M(I)*UN(SY)
2220 NEXT
2230 CLOSE 2:CLOSE 15:CB=1
2240 GOTO 370
2250 REM SAVE CURRENT SYSTEM TO DISK
2260 PRINT "[211]AVE CURRENT SYSTEM."
2270 INPUT"[212]YPE NAME OF FILE";A$
2280 IF LEN(A$)>13 THEN PRINT"[206]AME TOO LONG.":GOTO 2250
2290 OPEN 15,DV,15:C$=CHR$(13)
2300 OPEN 2,DV,2,"0:"+A$+".NB,S,W"
2310 GOSUB 3170
2320 IF ER<>0 THEN PRINT ER$(1);ER$(2);ER$(3);ER$(4):GOSUB 3160:GOTO 2230
2330 PRINT#2,SY;C$;NB
2340 FORI=1TONB
2350 PRINT#2,N$(I);C$;X0(I);C$;Y0(I);C$;Z0(I);C$;U0(I);C$;V0(I);C$;W0(I);C$;M(I)
2360 NEXT
2370 CLOSE 2:CLOSE 15:GOTO 370
2380 M(CB)=NV:GM(CB)=G*M(CB)*UN(SY):GOTO 370
2390 X0(CB)=NV:GOTO 370
2400 Y0(CB)=NV:GOTO 370
2410 Z0(CB)=NV:GOTO 370
2420 U0(CB)=NV:GOTO 370
2430 V0(CB)=NV:GOTO 370
2440 W0(CB)=NV:GOTO 370
2450 DT=NV:D2=DT*DT/2:D3=D2*DT/3:D4=D3*DT/4:GOTO 370
2460 REM DISPLAY CURRENT SYSTEM VALUES
2470 PRINTCHR$(147)"            [206]-[194][207][196][217] [211][201][205][213][204][193][212][207][210]"
2480 PRINTCHR$(176);L$;CHR$(174);
2490 L1$=CHR$(221)
2500 PRINTL1$" NAME: "N$(CB);TAB(79);L1$;
2510 PRINTL1$" BODY #"CB;TAB(17)"MASS:"M(CB);TAB(39)L1$;
2520 PRINTL1$" X:"X0(CB);TAB(60)"U:"U0(CB);TAB(79);L1$;
2530 PRINTL1$" Y:"Y0(CB);TAB(20)"V:"V0(CB);TAB(39);L1$;
2540 PRINTL1$" Z:"Z0(CB);TAB(60)"W:"W0(CB);TAB(79);L1$;
2550 PRINTCHR$(173);L$;CHR$(189)
2560 PRINT "NUMBER OF BODIES:"NB
2570 PRINT "TIME INTERVAL:"DT;TU$(SY)
2580 PRINT "UNIT SYSTEM: "UN$(SY)
2590 IF SP THEN PRINT "SPRITE MODE";CHR$(17):RETURN
2600 PRINT "HIRES POINT MODE";CHR$(17):RETURN
2610 :
2620 REM DISPLAY MENU
2630 PRINT "E[146]XIT       P[146]LOT       N[146]EW SYSTEM"
2640 PRINT "SC[146]ALE      D[146]ISPLAY    L[146]OAD [196][146]RIVE#"
2650 PRINT "S[146]AVE       F1[146] PREV BODY   F7[146] NEXT BODY"
2660 PRINT "X[146] POSITION Y[146] POSITION Z[146] POSITION"
2670 PRINT "U[146] VELOCITY V[146] VELOCITY W[146] VELOCITY"
2680 PRINT "M[146]ASS       T[146]IME       R[146]ENAME"
2690 RETURN
2700 REM PLOT POINT ON HIRES SCREEN
2710 IF SP=1 THEN 2750
2720 ML=HR+(YANDC8)*CF+(YANDC7)+(XANDC4)
2730 POKE ML,PEEK(ML)OREX(XANDC7)
2740 RETURN
2750 IF I>8 THEN RETURN
2760 X=X+24:Y=Y+50
2770 POKE VI+(I-1)*2,XANDC5
2780 POKE VI+I*2-1,Y
2790 IFX>C5THENPOKEHI,PEEK(HI)ORE2(I-1)
2800 IFX<256THENPOKEHI,PEEK(HI)AND(C5-E2(I-1))
2810 RETURN
2820 REM ML CODE FOR HIGH SPEED ERASE
2830 I=49152
2840 READMC:IFMC=256THENRETURN
2850 POKEI,MC:I=I+1:GOTO2840
2860 DATA173,52,3, 133,2, 173,53,3, 133,3
2870 DATA165,251, 160,0, 166,3
2880 DATA145,2, 236,55,3, 208,7, 166,2, 236,54,3, 240,9
2890 DATA230,2, 208,236, 230,3, 76,14,192, 96,  256
2900 PRINT CHR$(147);"[211]ELECT A SYSTEM OF UNITS."
2910 PRINT CHR$(17);"1. 1 PIXEL = 10^7 KILOMETERS"
2920 PRINT "   1 MASS = 1000 KILOGRAMS"
2930 PRINT "   1 TIME = 1 DAY"
2940 PRINT CHR$(17);"2. 1 PIXEL = 1 [193][213] (EARTH RADIUS)"
2950 PRINT "   1 MASS = 1 EARTH MASS"
2960 PRINT "   1 TIME = 1 DAY"
2970 PRINT CHR$(17);"3. 1 PIXEL = 1000 KILOMETERS"
2980 PRINT "   1 MASS = 1 KILOGRAM"
2990 PRINT "   1 TIME = 1 SECOND"
3000 PRINT: INPUT"[215]HICH SYSTEM";SY
3010 IF SY<1ORSY>3THEN370
3020 FORI=1TONB
3030 GM(I)=G*M(I)*UN(SY)
3040 NEXT
3050 GOTO 370
3060 :
3070 REM SWITCH PLOT SYSTEMS
3080 IF SP=1 THEN SP=0:GOTO 370
3090 SP=1
3100 FORI=15872TO15872+8*64:POKEI,.:NEXT:REM BLANK OUT SPRITE IMAGES
3110 FORI=0TO7:POKE15872+I*64,224:POKE15875+I*64,224:NEXT:REM FORM DOT SHAPE
3120 FORI=0TO7:POKE2040+I,248+I:NEXT:REM SET SPRITE POINTERS
3130 FORI=0TO7:POKEVI+39+I,I+1:NEXT:REM SET SPRITE COLORS
3140 POKE VI+29,0  :POKE VI+23,0  :REM COMPRESS SPRITES
3150 GOTO 370
3160 FORDE=1TO1500:NEXT:RETURN
3170 INPUT#15,ER$(1),ER$(2),ER$(3),ER$(4)
3180 ER=VAL(ER$(1)):RETURN
3190 CB=CB-1:IF CB<1 THEN CB=NB
3200 PRINTCHR$(19)CHR$(17)CHR$(17);B$;B$;B$;B$;B$:PRINTCHR$(19)
3210 GOSUB 2480
3220 FORI=1TO6:PRINTCHR$(17);:NEXTI:PRINTCHR$(29);:GOTO410
3230 CB=CB+1:IF CB>NB THEN CB=1
3240 GOTO 3200
3250 POKE53272,PEEK(53272)AND247
3260 POKE53265,PEEK(53265)AND223